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Abstract: Numerical solutions to fractional differential equations can be extremely compu- 



tationally intensive due to the effect of non-local derivatives in which all previous time points 
contribute to the current iteration. In finite difference methods this has been approximated using 
the 'short memory effect' where it is assumed that previous events prior to some certain time 
point are insignificant and thus not calculated. Here we present an adaptive time method for 
CN smooth functions that is computationally efficient and results in smaller errors during numerical 

simulations. Sampled points along the system's history at progressively longer intervals are 
assumed to reflect the values of neighboring time points. By including progressively fewer points 
as a function of time, a temporally 'weighted' history is computed that includes contributions 
from the entire past of the system, resulting in increased accuracy, but with fewer points actually 
^\ calculated, which ensures computational efficiency. 
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1. INTRODUCTION tribution of the memory effect associated with the his- 

(yQ tory of a system, which fractional methods must take 

(N Throughout the past few decades, fractional differential int o account. While this method can be applied to any 

^ equations have played an increasing role in the disciplines diffusion application modeled as a fractional process, it of- 

^ of physics (Compte (1997); Fellah et al. (2002); Metzler fers particular significance for modeling complex processes 

et al. (1999); Soczkiewicz et al. (2002)), chemistry (Goren- with long histories, such as lengthy molecular or cellular 

O flo et al. (2002); Seki et al. (2003)), and other engineering simulations. Approaches for numerically approximating 

O and science fields (Cushman and Ginn (2000); Baeumer solutions to fractional diffusion equations have been exten- 

^— 1 et al. (2001); Mathieu et al. (2003)). In particular, frac- sively studied (Chen et al. (2008); Langlands and Henry 

^ tional calculus methods offer unique approaches for mod- (2005); Liu et al. (2007); Mainardi et al. (2006); McLean 

• tH eling complex and dynamic processes at the cellular level. and Mustapha (2009); Podlubny et al. (2009); Yuste and 

Recent research using these methods have been applied to Acedo (2003)), but in general there is always a trade off 

H modeling ultrasonic wave propagation in cancellous bone between computational efficiency and the accuracy of the 

(Sebaa et al. (2006)), bio-electrode behavior at the cardiac resultant approximations. We developed a simple approach 

tissue interface (Magin and Ovadia (2008)), and describing that significantly reduces simulation times but which does 

spiny dendrite neurons through a fractional cable equa- not significantly affect the computed accuracy of the fi- 

tion (Henry et al. (2008)). Other work has shown that nal solution. This was accomplished through an 'adaptive 

endogenous lipid granules within yeast exhibit anomalous time step memory' method by changing the interval of 

subdiffusion during mitotic cell divison (Selhubcr-Unkel the backwards time summation in the Griinwald-Letnikov 

et al. (2009)). fractional derivative formulation used in many numerical 

schemes for discretizing fractional diffusion equations from 
Here we investigate the implementation and subsequent a cont i n uous Riemann-Liouville approach (Chen et al. 
computational performance of a fractional reaction-diffusion (2009); Yuste and Acedo (2003); Podlubny (1999)). In- 
equation (FRD) in numerical simulations. We introduce stead of incorporating every previous time point into the 
an 'adaptive time step' method for computing the con- most currcnt calculation, only certain points were recalled 

based upon their proximity to the currcnt point in time 

* This work was funded by NIH grant ROl NS054736. CM. was anc i weighted according to the sparsity of the time points. 

also partially supported through NIH T32 EB009380. 



This substantially reduced the computational overhead 
needed to recalculate the prior history at every time step, 
yet maintained a high level of accuracy compared to other 
approximation methods. We compared this new 'adap- 
tive time' approach with the 'short memory' principle 
described by Volterra (Volterra (1931)) and Podlubny et 
al (Podlubny (1999)) and examined differences in simu- 
lation times and errors under similar conditions for both 
methods. 

All calculations and simulations in the results were based 
on an explicit numerical method utilizing the Griinwald- 
Letnikov relationship for the fractional reaction diffusion 
equation given by: 



d 1 u{x, t) 



aV 2 u(x,t) + f(u) x = {xi, x 2 --.x N } (1) 



where 7 is the anomalous diffusion exponent, x is an N 
dimensional vector, a is the coefficient of diffusivity, and 
f(u) represents a consumption or generation term. 

2. METHODS 

2. 1 Griinwald-Letnikov formulation and simplification 

Consider the fractional diffusion equation with no reaction 
term, 



d~<u 
~dF 



= aS7 2 u 



(2) 



where a is the diffusion coefficient. Using the relation 
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the diffusion equation then becomes 
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~dt 
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where D 1 " 7 represents the continuous time Riemann- 
Liouville fractional derivative. We then implement the 
Griinwald-Letnikov definition for numerical simulations. 
This approach discretizes the Riemann-Liouville derivative 
through the relationship: 



From here, one can simplify the binomial coefficient to a 
simple recursive relation. 

Next, define a function ^(7, m) such that 



#y,m) = (-l)' 
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to 
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and substituting ( m 7 ) into eq. 6 yields 



(-l m )r(2- 7 ) 
m!r(2-7-m)' 



(9) 



Substituting the relation T(n) = (n — l)T(n — 1) into eq. 
8 results in 



■0(7,™) = 



(_l)m-l r(2 _ 7) 
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yielding the iterative relationship 



■0(7, m) = —'0(7, TO — 1) 



2 — 7 — to 
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For to = 0, 



(11) 
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This recursive function is valid for all 7 including sub- 
diffusion, standard diffusion and supcrdiffusion, so this 
equation is general over all regimes. Since there is no use of 
r(x), there are no values which tend to infinity or limits to 
compute. In addition, since these values are used over and 
over again during the course of a numerical simulation, 
they can be precomputed for values of to = to to = N, 
where N is the number of time points, resulting in a signif- 
icant savings in performance. A similar simplification has 
been previously discussed and used in Podlubny (1999). 

2.2 Discretization 

The substitution of tp(j,m) into eq. 7 results in 



D»u = ^£{-ir(£\f{t-mT) (5) 

m=0 ^ ' 

where (^) represents the binomial coefficient given by Then, using 

n\ r(n + l) 



m!(n — to)! m\T(n + 1 — m) 
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Applying eq. 5 to the fractional diffusion equation (eq. 4) 
yields 
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as the formulation in two dimensions, one can discretize 
the function into a finite difference based FTCS scheme 
(forward time centered space) on a grid l where k = 
t/A t ,j = xi/A x ,l = x 2 /A X} using the relations (Press 
et al. (1992)) 
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In the discrete limit where t — >• A* 



u j,i u j,i a; 



(16) 



x m=0 



where 8 k L m is the finite difference kernel given by 



(17) 

Adding a consumption/generation term is straightforward 
in this implementation. For example, take an exponential 
decay term given by 



du 
~dt 



(18) 



with the complementary finite difference relation 
v k+1 - v k 

■ — = -0u k iA t 



Incorporating eq. 18 into the form of eq. 4 results in 
du cP" 1 



(19) 



dt dn 
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which gives the full finite difference implementation in two 
dimensions 
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3. RESULTS 

3.1 Diffusion simulation examples 

Fig. 1 shows the results of simulating the equations for 
four different values of 7. Simluations were run on a 100 
x 100 grid with initial conditions L/50,50 = 0.1, 50 = 



77° - r/° 

u 50.51 — u 4< 



49.50 



l/K, 



0.05 and zero elsewhere. 



Dirichlct boundary conditions were implemented and the 
simulation edge set to zero. In these simulations, a = 1, 
13 = 0, A t = 0.5, A x = 5, and run until t = 100. 

3.2 Adaptive Memory Method 

In eq. 21 each iteration requires the re-calculation and 
summation of every previous time point convolved with 
V>(7,m). This becomes increasingly cumbersome for large 
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Fig. 1. Simulation results for 7 = 0.5, 0.75, 0.9, 1.0 (for 
traces from top to bottom) in one dimensional space 
(panel A) and time (panel B). As the subdiffusion 
regime is approached the profile becomes more and 
more hypergaussian. 

times, which require significant numbers of computations 
and memory storage requirements. To address this, Pod- 
lubny et. al (Podlubny (1999)) introduced the 'short mem- 
ory' principle which assumes that for a large t the role 
of the previous steps or 'history' of the equation become 
less and less important as the convolved ip(j,m) shrinks 
towards zero. This would then result in approximating eq. 
21 by truncating the backwards summation, only taking 
into account times on the interval [t — L, t] instead of 
[0, t], where L is defined as the 'memory length' (eq. 7.4 in 
Podlubny (1999); Fig 2). While computationally efficient, 
this approach leads to errors in the final solution since 
not all points are counted in the summation. Despite 
the resultant errors, this numerical method represents a 
powerful and simple approach for providing a reasonable 
trade off between computational overhead and numerical 
accuracy. In the context of the implementation derived 
here, it would result in a discretization scheme given by 
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As an alternative to the method of Podlubny, we propose 
an adaptive time approach that also shortens computa- 
tional times but which at the same time results in much 
greater accuracy than the use of 'short memory'. We 
achieve this by introducing the concept of an 'adaptive 
memory' into the Griinwald-Letnikov discretization. 

The underlying principle of the adaptive time approach 
is that relative to the current time point previous time 
points contribute different amounts to the summation. 
Values relatively closer to the current time point will have 
a greater contribution to the current numerical calculation 
than values many time points back due to the multiplier 
tpil, m). For smooth functions, as m increases and ipij, m) 
decreases, neighboring points in the summation exhibit 
only small differences. Consequently, one can take advan- 
tage of this and utilize an 'adaptive memory' approach in 
which neighboring values at prior time points are grouped 
together for defined increments and the median taken 
as a representative contribution for the increment and 
weighted according to the increment length to account 
for the skipped points. This results in fewer time points 
and fewer computations in the summation. Algorithmi- 
cally, for an arbitrary a time steps back from the current 



time point k for which the history of the system is being 
computed, consider an initial interval [0, a] for which all 
time points within this interval are used in the summation 
and therefore contribute to the the Griinwald-Letnikov dis- 
cretization. Let subsequent time intervals become longer 
the further away from k they are and within them only 
every other d time points are used in the summation 
term, i.e. only the median values of a temporally shifting 
calculation increment of length d along the current interval 
are considered. As subsequent intervals become longer so 
does d, thereby requiring less points to compute. 

Definition 3.1 Let k be the current iterative time step in a 
reaction diffusion process for which a Griinwald-Letnikov 
discretization is being computed. Consider an arbitrary 
time point a in the history of the system backwards from 
k. For i e Ni, i ^ 1, define an interval of this past history 

by 

/ = [a 1 - 1 + i, a 1 } (23) 

where Ni represents the set of natural numbers beginning 
from one. Given how the indices i are defined, the very first 
interval backwards from k is independent of 23 and is given 
by [0, a]. This is considered the base interval. Subsequent 
intervals arc defined as a function of this base, i.e. as a 
function of a and determined by 23. 
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Let imax be the value of i such that k £ I ma x — [a lma!r 1 + 
'maua 1 """]. The complete set of intervals then is defined 
as ( = {/ = [a 1 - 1 +i,a l ] : i eNi,i ^ l,i < i max }. 

Definition 3.2 For the set of intervals C as defined in 
Definition 3.1, D = {d = 1% - 1 : i G Ni,i ^ 1} is the set 
of distances d by which the corresponding intervals in Q 
are sampled at. 

Proposition 3.1 In general, for two dimensional diffusion 
without consumption or generation terms for any interval 
as defined in Definition 3.1, the Griinwald-Letnikov dis- 
cretization with adaptive time step memory is given by 
l 

^(7,»)$r+--- ( 24 ) 
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where p € Ni, and for each i (i.e. for each interval) M = 
{irii = a 1 +(2i— 1)77 — i+1 : 77 € Ni & m., < m max } is the 
set of time points over which tp("/,m)SjJ' m is evaluated. 
Since the time point k may be less than the full length 



m, 



< 



\m k - 



of the last interval I max 
maximum value in I max that is evaluated 
element in the set M for I m ax- 



represents the 
i.e. the last 



Proof The first summation represents the basis interval 
and the limits of the summation operator imply the 
contribution of every time point, i.e. n € Ni. For intervals 
beyond a: any arithmetic sequence defined by a recursive 
process u v = + d, rj € Ni for some distance d, the 
if h value in the sequence can be explicitly calculated as 
^77 = v\ + (jj — l)d given knowledge of the sequence's 
starting value v\. For the set £ this implies that v\ = a' _1 + 
i and d — 2i — 1 for a given i. This then yields v v = a ,_1 + 
i + (77 - l)(2i - 1) = a^ 1 + (2i - l)r) - i + 1 := m ; as 
required. The outer summation collects the summations of 



Fig. 2. Short memory and adaptive memory methods for 
estimating the the Griinwald-Letnikov discretization. 
Both approximations rely on the sharply decreas- 
ing function ^(7, m) as m increases to omit points 
from the backwards summation. While short memory 
defines a sharp cut off of points, adaptive memory 
provides a weighted sampling of points for the entire 
history of the system. Points included in the compu- 
tation by each method are highlighted in red. See text 
for details. 

all intervals that belong to the set £■ The last summation 
on the interval [m max +i, k] ensures that the final point (s) 
-Df the backwards summation are still incorporated into the 
p computation even if the arbitrarily chosen value of a does 
not generate a final increment length that ends on k. □ 

Note that D is not explicitly needed for computing equa- 
tion 24 because the distances d are implicitly taken into 
account by £. Using the median value of each increment 
avoids the need for interpolation between time points. The 
implementation of the adaptive memory method described 
here is necessarily limited to smooth functions due to the 
assumption that the neighbors of the median values used 
in the summation do not vary much over the increment 
being considered (but see Discussion below) . This method 
essentially allows for a contribution by all previous time 
points to the current calculation, yet reduces computa- 
tional times by minimizing the total number of calcula- 
tions. For the interested reader the full source code for 
the implementation is available online from the authors' 
website at www.silva.ucsd.edu. 

The results of using various L (for short memory) and 
interval steps (for adaptive memory) are shown in Fig. 3. 
Increasing the values of L and a resulted in a decrease 
in the error of the estimated results but at the cost of 
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Fig. 3. Comparison of the error between adaptive memory 
and the short memory as a function of the calculation 
time (x-axis: computation run time in seconds) ex- 
pressed as a percentage error relative to the computed 
value for the full non-shortened approach (y-axis). 
Four different values of 7 are shown. 

increased computation times. Conversely, decreasing the 
values of L and a resulted in a decrease in computation 
times, but at the expense of accuracy. In all cases however, 
the adaptive time method had a significantly lower error 
for the same simulation time, and also reached a clear 
asymptotic minimum error much faster then the minimum 
error achieved by the short memory method. In these 
simulations, a = 1, /3 = 0, A t = 1, A x = 10, using a 
20 x 20 grid, and ran until t = 1500 where Uf 010 = 10. 

The error for the 'short memory' method increased com- 
paratively quickly and got worse as 7 — > 1. This was due 
to the fact that the diffusion from the initial point source 
was initially quite fast near t = 0, which were the first 
time points cut by the 'short memory' algorithm. Because 
the component near t — represents the majority of the 
immediate past history function, large errors result even 
when L — > t. 



4. DISCUSSION 

The regular diffusion equation has been a cornerstone of 
transport modeling for decades. Regular diffusion relies 
on the assumption of a Markovian process that is not 
necessarily present in natural systems. 

One approach to modeling these non-markovian processes 
is using the fractional diffusion equation (cq. 1). Math- 
ematically such methods have been around for decades 
but it is relatively recently that they have been applied 
to the natural sciences. Computationally efficient methods 
for numerically evaluating these equations are a necessity 
if fractional calculus models are to be applied to modeling 
real world physical and biological processes. 



It should be noted that while in this work the simulations 
were done in the subdiffusive regime for a simple point 
source, the methods we derive here are directly applica- 
ble to more complex sources or superdiffusion (7 > 1). 
However, complex fast-evolving past histories in the form 
of a forcing function (/(it)) or oscillations generated in a 
superdiffusion regime will result in much larger errors for 
both the short and adaptive memory methods. In the case 
of the adaptive memory method introduced in this paper 
this is due to its 'open-loop'-like algorithm that blindly 
increases the spacing between points in the summation 
as the calculation goes further back in time and which 
does not take into account the speed of evolution of the 
equation. Adaptive time approaches for regular differential 
equations often make the integration step size a function 
of the derivative, i.e. more closely spaced time steps are 
used when the function is oscillating quickly and cannot 
be averaged without loss of accuracy, and widely spaced 
time steps are used where the function is smooth. In the 
current implementation we have assumed that the past 
history function ^(7, m)S'lJ m is smooth. Subsequent work 
will implement a 'smart' 'closed-loop'-like algorithm where 
the step size of the backwards summation is dependent on 
the derivative of the past history function, i.e. a form of 
feedback. This would optimize the speed of the simulation 
while reducing the error due to the averaging of time 
points in the backwards summation, ultimately resulting 
in low errors for both high frequency forcing functions in 
the subdiffusion regime and for oscillations intrinsic to the 
superdiffusion regime. 
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